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ABSTRACT 

Context. The smoothed particle hydrodynamics (SPH) technique is a well-known numerical method that has been applied to simulate 
the evolution of a wide variety of systems. Modern astrophysical applications of the method rely on the Lagrangian formulation of 
fluid Euler equations which is fully conservative. A different scheme, based on a matrix approach to the SPH equations is currently 
being used in computational fluid dynamics (CDF). These matrix formulations achieve better interpolations of the physical magnitudes 
but they are, in general, not fully conservative. The matrix approach to the Euler equations has never been used in astrophysics. 
Aims. We develop and test a fully conservative SPH scheme based on a tensor formulation that can be applied to simulate astrophysical 
systems. 

Methods. In the proposed scheme, derivatives are calculated from an integral expression that leads to a tensor (instead of a vectorial) 
estimation of gradients and reduces to the standard formulation in the continuum limit. The new formulation improves the interpolation 
of physical magnitudes, leading to a set of conservative equations that resembles those of standard SPH. The resulting scheme is 
verified using a variety of well-known tests, all of them simulated in two dimensions. We also discuss an application of the proposed 
tensor method to astrophysics by simulating the stability of a Sun-like polytrope calculated in three dimensions. 
Results. The proposed scheme is able to improve the results of standard SPH in the two-dimensional tests, especially in the simulation 
of subsonic hydrodynamic instabilities. Our results for the stability of the Sun-like polytrope suggest that the new method can be used 
in astrophysics to carry out three-dimensional calculations with a computational cost that is only slightly higher (i.e. < 50% for a 
serial code) than that of a standard SPH formulation. 

Conclusions. A formalism based on a matrix approach to Euler SPH equations was developed and checked. The new scheme is 
more accurate because of the re-normalization imposed on the interpolations, which is fully conservative and probably less prone to 
undergo the tensile instability. The analysis of several test cases suggest that the method may improve the simulation of both subsonic 
and supersonic systems. An application of the tensor method to astrophysics is, for the first time, successfully carried out. These 
encouraging results indicates that more work should be invested in the applications of matrix SPH formulations to astrophysics. 
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1. Introduction of all, the numerical noise is usually larger than in other tech- 
niques, which poses a problem for applications involving very 

The hydrodynamic method known as smoothed particle hydro- subsonic moV ements, as in the case of hydrodynamic al instabil- 

dynamics (SPH) is a grid-less Lagrangian approach tocontin- ities Another difficulty has to do with the use of the artificial 

! ™™ mechanics devised by |Gingold & Monaghan| OSZZb and viscos i ty (A V) formalism to handle shock waves and smooth the 

iLucyJ (11921. As demonstrated in an endless amount of apph- inhere nt presence of noise. Several solutions and recipes have 

■ cations to computational fluid dynamics, ranging from large- been invoke d to handle these problems, none of t hem totally sat - 

scale astrophysical simulations to nuclear physics on the Fermi is f actorVj among them: variable AV coefficients dMorris II 19971). 

level, SPH is a robust and confident way to simulate the dy- ^fc^i hea t fluxes to smooth pressure discontinuities dPricTl 

namical evolution of deformed systems. The success of SPH 2 008), non-standard approximations to the momentum equation 

relies strongly on the way gradients are estimated. For exam- jAbdlEoTlh . and a refined treatment of kernel interpolations as 

pie, the value of density in a given spatial coordinate is obtained in the tec hnique known as moving i east squares interpolation 

from particles located in the neighborhood using an interpolat- (j e jviLS" iDilts l!T999h 
ing function called kernel. As the kernel is an analytical differ- 

entiable function, it is easy to obtain gradients just evaluating l n this paper we present an approach to the momentum and 

the gradient of this weighting kernel function. Then it is quite energy equations that combines a novel way to estimating gradi- 

straightforward to write the Euler equa tions of fluid mechanic s ents and the variational principle, and allows an improved mod- 

in terms of the kernel and its derivatives (iMonaghanl 19921 12005). eling of physical phenomena, especially in systems subjected to 

Despite the success, SPH has also a number of shortcom- small perturbations. The proposed formalism, based on an in- 

ings and weak points that make the technique less useful than it tegral approach to the derivatives (IAD), includes the standard 

should be to handle a number of well-identified problems. First formulation as a particular case. It leads to a matrix formulation 
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similar to that of MLS methods, but the form of the SPH equa- 
tions in IAD remains closer to that of standard SPH. The text 
is organized as follows. In Sect. [2] the mathematical formalism 
for calculating gradients is described and discussed. In Sect. [3] 
we develop a formulation of the momentum equation compat- 
ible with the variational Euler-Lagrange principle in the SPH 
framework, and propose a consistent and fully conservative en- 
ergy equation. Section [4] is devoted to validating the hydrocode 
through the simulation of four known tests where standard SPH 
has, to a smaller or larger degree, difficulties. In Sect.[5]we elab- 
orate an extension of the method that by combining the exact 
estimation of the derivative of linear functions and good mo- 
mentum and energy conservation, can be applied to simulate lin- 
earized systems. The ability of the proposed method to simulate 
the structure of real three-dimensional (3D) astronomical bod- 
ies is discussed in Sect. [6] Finally, the main conclusions of our 
work, as well as some comments about the shortcomings of the 
developed scheme and future lines of improvement are outlined 
in the last section, which is devoted to conclusions. 



2. Integral approach to first derivatives 

There is an ample literature desc ribing the treatment of the 
first a nd second derivatives in SPH (Monaghan 2005; Rosswog 
l2009h . In the following we take a different route to estimating 
first derivatives which leads to a more general expression than 
the most common extant derivative procedures. Our method re- 
sembles to those known as MLS methods where interpolations 
relative to gradient estimation were constrained to enhance lin- 
earity by imposing a minimization procedure on th e errors. 

Inspired by the treatment of second derivatives (Brookshaw 

Il985h we define an integral expression from which we will es- 
timate gradients. As in the case of second derivatives, it is ex- 
pected that an integral approach will lead to a less noisy estima- 
tion of first derivatives. Therefore, we define 



/(r)= f [/(r')-/(r)](r'-r)W(|r'-r|,^r' 3 , 
Jv 



(1) 



where /(r) is any differentiable function and W(\r' - r\,h) is a 
spherically symmetric interpolating kernel, usually a sharp-like 
Gaussian of width h (where h is usually referred as the smooth- 
ing length). The size of h is often interpreted as the local reso- 
lution of the simulation. Expanding the factor [/(r') - f(vj\ to 
first order 



/(r') - /(r) = Vf • (r' - r) + higher order terms (2) 
and writing 



h = j [/(r') - /(r)] (x' k - x k )W(\r' - r|, h)dr' 3 ;k=l,3. (6) 

Once the linear system in Eq. © is solved we get the value 
of the gradient of the function /. In the case of spherically sym- 
metric kernels, t\ \ — T22 — T33 and T{j - t y,- = 0, i + }■ Taking 
for instance the Gaussian kernel 



W G (u, h) = — exp(-M 2 ) , 



(7) 



where n is the dimensionality of space and u = ^ r h r ^ , it leads 

to Tii = y • Therefore, for the Gaussian kernel the gradient of / 
reduces to 



Vf : 



)-/(r)] VWflr' -r\,h)dr' 3 



(8) 



which is the same expression as Eq. (2. 13) of Monaghan (2005) 
even though it has been obtained through a different procedure. 
Since the standard kernels are spherically symmetric Eq. (0 can 
be simplified to 



Vf = J /(r')VW(|r' - r|, h)dr' 3 , 
which in SPH parlance becomes 

V a f = Y — f b VV? ab (h a ), 

V Pb 



(9) 



(10) 



which because of its simplicity and good performance has be- 
come the most popular procedure for calculating first derivatives 
in SPH. Therefore, Eq. ( fTOb is a particular case of IAD, Eq. (0). 

Nevertheless, many of the above symmetry properties of ten- 
sor T = {Tjj} are lost when integrals are converted into finite 
summations. For example, the elements on the diagonal of ma- 
trix T in Eq. <j4j do not necessarily have the same value and 
elements outside the diagonal can differ from zero. However, 
the formulation of the first derivative in terms of the matrix Eq. 
© has a very interesting feature: the derivative of a linear func- 
tion is always exact by construction. The demonstration of such 
property is straightforward in ID SPH: 



df 
dx 



Zb'^(fb-fa)(x h -X a )W ab 



(11) 



taking f a = px a + q,fb = px b + q the expression above gives 

{%)a = P- 

In two dimensions, the explicit solution for Eq. is written 



as 



(r' - r) = {x\ - x{)i + (x' 2 - x 2 )j + (x' 3 - x 3 )k , 
equation ([!} can be expressed as a matrix equation: 



where 



and 



df/dxi 




tu 


T\2 


Tl3 ' 


— j 


\h 1 


df/dx 2 




T21 


T22 


T23 
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df/dx 3 




. T31 


T32 


T33 . 







= f(xl - xdix'j - Xj )W(\r' - r|, h)dr' 3 ; i, j = 1 , 3 



(3) 



(4) 



(5) 



df_ 
dx\ 



■ cnh + c n h ; 



where 



dx 2 



C2\h + C22I2 : 



(12) 



22 = (t 22 - ^) • 



(13) 



Unfortunately, direct application of the tensor scheme above 
leads to SPH equations of movement that are not fully conserva- 
tive. Nevertheless, a slight modification of the scheme makes it 
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possible to build exact conservative equations by taking advan- 
tage of the symmetry properties of the kernel to discretize the 
vector I(r) as 



/« = j [f(r') - /(r)] (r' - r)W(|r' - r|, h)dr" - 



X — /(r b )(r b - r a )W(|r b - r a |, h a ) . 



(14) 



Thereafter, the formulation that results after substituting Eq. 
( TPfl i into Eq. (|4]i is labeled IADq. This restricted interpretation 
of Eq. (H]i implies that exact evaluation of the derivative of lin- 
ear functions cannot be achieved if we are to also preserve the 
exact conservation of momentum and energy. Nevertheless, we 
present below strong indications that using IADq leads to a better 
evaluation of the derivative of linear functions than the standard 
procedure given by Eq. ( TTOb . especially when a small or mod- 
erate number of neighbors is used to compute summations. We 
can illustrate this point in the following numerical experiments. 

As a first test, we simulate a static system where we compare 
the gradient of a linear distribution of density in one dimension 
obtained with IAD, IADq and standard schemes as a function of 
the number of neighbors rib- A linear density profile p(x) = 1 + 
x was obtained aligning N= 100 equidistant particles of adequate 
mass along the x-axis. The density was calculated using 



Pa = ^mWahQr), ■ 



A). 



(15) 



Note that even though IAD should provide the precise 
derivative, in practice it never does owing to the small errors in 
the calculation of density. It is also worth noting that the standard 
calculation of the derivative has two additional potential sources 
of error. For a linear density profile pb - p a + p(xb - x a ), and 
using Eq.dTOb. we have 



Id \ STD 

/ = 7! —paSaWab + p V j — (x b - X a )V a W ab 

\ dx L V Pb 



(16) 



= e\ + e 2 p , 



where p is the 'exact' derivative and e\ 0, ei — 1. On the other 
hand, there is only one source of error, 63 - 0, if IADq is used, 
given by 



IAD ° 'fPaiXb ~ X a )W ah 



Xbf^Xb-XafWab 



+ p = e 3 + p . 



(17) 



In Fig. Q] (upper and bottom left) we represent the value of 
the relative error in the derivative e = |^ - p\/p, with respect 
to the analytical value p as a function of the smoothing length h, 
normalized to the inter-particle distance A, for the Gaussian and 
cubic spline kernels respectively. Independently of the kernel, 
the error when full IAD linear interpolation was used is always 
much smaller than for the other methods. When IADq or the 
standard derivative were used, the error increased appreciably. 
The error is large when the smoothing length is shorter than the 
inter-particle distance but decreases rapidly when h increases, as 
expected. Nevertheless, the error in the standard calculation al- 
ways remained larger, especially for a small or moderate number 
of neighbors («/, < 5) regardless of the kernel we used, but es- 
pecially for the cubic spline. Since the particle sample is highly 



ordered, the error curve of the Gaussian kernel is much smoother 
than that of the cubic spline owing to the infinite range of the ker- 
nel. We verified that when the Gaussian is truncated to 2h, the 
error profile becomes qualitatively similar to that of the cubic 
spline. However, the error profile for IADq follows almost ex- 
actly the error curve for the gradient of density calculated as a 
simple quotient for adjacent particles. In this sense, IADq 
seems to be more coherent with the computed density distribu- 
tion of the sample than IAD, a trend that also holds in 2D, as 
shown below. 

The contributions to the total error of €\ , €2 and 63 in Eqs. ( fT6b 
and ( TTTl i are shown in the rightmost part of Fig. 1. The contribu- 
tion of e\ and £j are similar but the error €2, which affects only 
the standard scheme and exhibits a large oscillation in the case 
of compact supported interpolators, represents the main channel 
to the total error in the derivative. In general, tensor methods 
become more reliable for higher dimensionality. Thus, we con- 
ducted a similar numerical experiment in more than one dimen- 
sion where the matrix approach can potentially be more benefi- 
cial. 

As a second test, a bi-dimensional system was set us- 
ing a sample of N=62500 equally spaced particles on a two- 
dimensional lattice. The mass of the particles along the jt-axis 
was conveniently modified to reproduce a linear density profile, 
p(x, y) = 1 + x. The first derivative of the density,^ for different 



number of neighbors was obtained using Eqs. (|4]l and (TToT >. In 
the upper row of Fig. |2]we show the relative error in the deriva- 
tive, e = \£ - p\lp, with respect to the analytical value p as a 
function of the smoothing length. Calculations were carried out 
using the Gaussian and the cubic spline kernels. The trend here 
is similar to that of the ID calculation but magnitudes fluctuate 
less, especially in the case of the cubic spline. When the number 
of neighbors is large the error is small in all three methods, but 
tensor schemes are more accurate than the standard procedure. 
Nevertheless, for shortening smoothing lengths, h < 1.5A, the 
error curve steepens and for h < A standard SPH begins to give 
unacceptable results while the tensor estimations of the deriva- 
tive remain applicable. This result also agrees with the ID calcu- 
lation shown in Fig. 1, which shows that working with «/, < 3,4 
in ID may give an appreciable error in the derivative of a lin- 
ear function if Eq. ( TTOb is used to compute the derivatives. As 
in the ID case, the gradient calculated as a direct fraction Ei-^s. 

& x b~ x a 

(crosses in Fig. 2) remained closer to the results obtained using 

IADq. 

Finally, we considered a spherically symmetric linear pro- 
file p(x, y) — 1 + r, where r = ^jx 2 + y 1 . In the bottom row of 
Fig. [2] the value of the relative error e is shown as a function of 
h(A). The profile of the error follows a slightly different trend 
than in the previous case. As before, the error is large when the 
smoothing length is shorter than the inter-particle distance but 
in the case of the Gaussian kernel there is a minimum in the er- 
ror profile at h 1.1 A. From that point on the error smoothly 
increases. This behavior is not so clear in the case of the cubic 
spline where the error curve roughly stabilizes above h = 2A. 
There is, therefore, an optimal value of h in this case that min- 
imizes the errors in the first derivative. Close to the center lin- 
earity is lost because of the symmetry imposed on the distribu- 
tion. When the value of h is large, interpolations are affected by 
that geometrical constraint and, at some point, the advantages of 
working using a large smoothing length are lost. Nevertheless, 
the errors calculated using JADo are always smaller than that of 
the standard scheme. Again the gradient of density estimated as 
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a simply quotient yzp is more closely reproduced by IADq than these elements in cartesian 2D were given by Eq. ( fT3l l. Equation 
by the other schemes. (l24l can be used to directly compute the second term on the right 

of Eq. ([20]). 

We can evaluate, in the same way, for particle b 

3. The momentum and energy SPH equations using 

IAD o dp b 4^ 

o — = / , / , m c c ijtb (h b )(x u - xj,b) W bc (h b ) . (25) 

Nowadays the most accurate formulation of the momentum <JX^ b j 1 *^ 
equation in SPH comes f rom the variational principle. It has been 

shown ( Monaghan 2005, and references therein) that the solution To calculate Jp- we use Eq. ( fl"5l ): 
of the Euler-Lagrange equations leads naturally to a very sym- 
metric scheme including the effects of spatial gradients on the dp b _ dW ab (h b ) 
smoothing length. The resulting equation conserves momentum dx, a " dxj b 
by construction and, equally important, ensures perfect energy . . 

conservation for non-dissipative systems. It is shown below how The last term ln E 1- ® als ° a PP ears d ™ n S the evaluation of 

the tensor approach built using the IADq scheme can also be ^bPb, namely 
compatible with the variational principle. 

The Lagrangian of the system is evaluated as d dW bc (h b ) dW ab (h b ) 

"3 — = A me — 3 = ... + m a — + ... (27) 

'I \ oxij, ^ ox Ub dxu, 



^ m \^v 2 b - u b (p b , s b )j , (18) 



Comparing Eqs. d27l i. (l26l i and d25l l. we finally get 



where v b , u b and s b are the velocity, specific internal energy and 

specific entropy of particle b. Using this Lagrangian and assum- dp b 

ing isentropic evolution, the Euler-Lagrange equations subjected q x . ~ a b ^ J< b x h a > ab( b) • 



to the constraint p/i 3 = constant, lead to the equation of move 

Pb 



7=1 



ment (ISpringel & Hernquist ll2002H T l Substituting Eqs. (|24) and (g8j into Eq. ©, the /-component of 

the momentum equation for particle a is given by 



Zf b 
m b —r—^aPb, (19) 

where Q b = (1 - dp/dh £ c m c dW bc /dh) is a term accounting *;, a = - ^ m b\ — - 2 ^i,ah(h a ) + — — jM' iab {h b ) 

for the gradient of h. The /-component in Eq. (JT9J can be written b=i ^ "P" il bP b 

V ff, dp,, P fl <9p fl where 

w a x ljfl = - > m fe -r— m a -r— - — . (20) 

pfah dx i>a piQ. a dx Ua 

An estimation of the density gradients using the tensor Eq. © is . v -1 w -.n, 

&i,abQ l <d = 2^ c uA h a)( x j,b- Xj, a )W ab (h a ) , 



(29) 



(30) 



Vp = CI, (21) 



7=1 

d 



where C = T . Elements of matrix 7" are those defined in Eq. W iab (h b ) = Y c u , b (h b )(x lb - x M )W ab (h b ) . (31) 

(O after changing integrals to summations. For particle a, they 
are defined to be 

mh It is then straightforward to derive the appropriate energy equa- 

Tij, a = /. — {*i,b ~ Xi,a)(*j,b - Xi a ) W a b(K) , (22) tion 

and the y'-component of vector I in the IADq approach is . . n b d 

Z ~TT = ^ " 7 A A m 6(v/, fl - v ifi )^ ab (h a ) . (32) 
Xj^WWAo). (23) \dt) a Cl a plj^ x j^ 

b 

Thus, for particle a the /-component of density gradient is Since the magnitudes ^ ab , defined by Eqs © and CD, are 

antisymmetric, the conservation properties of the IADq formula- 
tion are identical to those of standard SPH. 



Qp a ^ As in standard SPH, an artificial viscosity term has to be 

-q^ — = 2_j c 'j,a(h a ) Ij,a(ha) - added to stabilize the numerical scheme and handle shocks. The 

''" i=i AV increases the effective pressure in those zones of the fluid 

d n b ' that undergo compression. Including the viscous acceleration, 

2 J] m b c ila (h a )(Xj >b - Xj, a ) W ah (h a ) , the momentum equation is given by 



j=\ b 



where c, ; fl are the elements of matrix C associated with particle ^ ( p a p h 

a, and d is the dimension of the space. Explicit expressions for x 'a = — / . n% — j ^i,ab{h a ) + — j^Lab^b) + H ab 

' b= i V'^aPa ^ l bP b 

In that paper the Vh contribution is estimated using the Lagrange (33) 

multipliers technique. where Yl ab accounts for the viscous pressure 
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n,,/, - < — S for Vab ■ y " h < °' (34) 

otherwise, 

where the symbols have their usual meaning and r a b means 
r fl - r^,. The coefficient fi a b is 



Hab - — 



hobnob ■ Vab 



(35) 



To preserve momentum conservation, we evaluate the arithmetic 
mean of J[ 



&i,ab - ^ [&i,ab(ha) + & Lab(hb)] 



(36) 



The inclusion of artificial viscosity in the energy equation 
leads to 



f du\ 



n b d 

V V 



dt 



la b=\ i=l \ Ll aPa L 



ab 



(37) 

As a representative example of the ability of the new scheme 
to handle subsonic physics, we consider the stability of an in- 
homogeneous system in mechanical equilibrium. A sample of 
N=62 500 particles were evenly distributed within a square lat- 
tice with periodic boundary conditions ensuring that the initial 
density was set to po = 1 g.cm~ 3 . The density was then per- 
turbed at random allowing maximum percentage variations of 
5% across the lattice. The internal energy of each particle was 
adjusted to ensure isobaricity with Pq — 1 dyn.cm 2 leading to a 
sound speed value c s =* 1.3 cm.s~'. The evolution of the system 
was followed during almost a sound crossing time, t s - 0.5 s 
using the three schemes, IAD, IADq and standard SPH, with 
rib = 30 and 100. Calculations using full IAD were carried out 
using the Eqs. ( 1521 and ( f53b described in Sect. [5] Hereafter, the 
word standard, STD, refers to the modern, fully conservative, 
Lagrangian formulation of SPH that explicitly includes the gra- 
dient of t he smoothing le ngth parameter in the scheme (see for 
instance iRosswog I (120091) and references therein). The results 
of the simulations are summarized in Fig. [3] where we present 
the evolution of the standard deviation in the pressure P(t) with 
respect to Po for the three aforementioned schemes. For this par- 
ticular test it is clear that both tensor schemes provide a more 
persistent mechanical equilibrium than standard SPH during a 
sound crossing time, especially when full IAD is used. As ex- 
pected, the standard deviation in the pressure decreases as the 
number of neighbors increases. 



3.1. Computational issues 

The increase of the computational requirements of the integral 
method is more evident now. With respect to memory require- 
ments, the nine coefficients (in 3D) of matrix 7~ and those be- 
longing to matrix C can share the same memory space. Thus, 
nine coefficients per particle have to be stored to compute Zft^ob 
in Eq. (f33). While the increase in memory requirements are not 
a significant problem, the burden in the computational time is 
larger. First of all, one has to compute six of the nine coefficients 
of symmetric matrix T, Eq. (T22l . Although these calculations in- 
volve simple operations, they have to be performed in a separate 
tree-walk because previous knowledge of density is necessary. 
Afterwards, matrix T has to be inverted in order to calculate the 



coefficients cy jfl for particle a. Fortunately, these matrix inver- 
sions do not take much time because they can be completed out 
of the tree structure. Finally, one has to compute the nine quan- 
tities of coefficients tti,ab associated with each pair of neighbors 
particles a, b; instead of the three that are necessary in the stan- 
dard formulation, and use them to compute both the momentum 
and the energy equations. 

It is obvious that the proposed method demands a larger 
computational effort than the standard SPH. Nevertheless, the 
extra burden could be small if the physical problem under simu- 
lation requires the calculation of long-range forces (i.e. gravity) 
or involves complex physics (i.e chemical or nuclear reactions, 
transport phenomena). In addition note that, according to Figs. 
Q~|and[2] for standard SPH to achieve a similar accuracy as IADq 
in calculating the first derivative, the number of neighbors has 
to be larger. To maintain the spatial resolution, the total number 
of particles also has to be larger by the same factor. Therefore, 
changing both the number of neighbors and total number of par- 
ticles could have a larger impact on the computational perfor- 
mance than using IADq. 

4. Hydrodynamic tests 

The aim of this section is to check the theory described in Sects. 
|2]and[3]using four tests, all of them widely used in CFD. The first 
two simulations concern the growth of the KH and RT hydro- 
dynamic instabilities, whose evolution always remains subsonic 
and, as we demonstrate, the fully conservative IADq scheme 
leads to improved results with respect to the standard method 
for the same identical initial conditions. The last two tests are 
related to strong supersonic flows where shock waves take over, 
as in the wall heating shock and Sedov tests, where the use of the 
tensor route to estimate derivatives again provides more accurate 
results, especially for the wall heating shock test. 

The simulations described in this section were carried out 
in a cartesian two dimensional scenario because of the higher 
achieved resolution, easiest initial setting and analysis of the re- 
sults. Periodic boundary conditions were imposed in all tests 
except in the RT case, where particles from then top and the 
bottom of the box were fixed. A perfect gas equation of state 
(EOS) P — (y- l)pu, where u is the specific internal energy and 
y = 5/3, is always used. The coefficients a, ft of artificial viscos- 
ity, in Eq. d34l , were set to 1 and 2 respectively, except for the 
wall heating test where we took a = 1.5, /? = 3. The cubic spline 
was used to perform the SPH interpolations but several calcula- 
tions carried out with other interpolato rs (i.e. harmonic kernels 
with index n=3: ICabezon et all ([2008)) did not lead to signifi- 
cant differences. We used a self-adaptive smoothing length pa- 
rameter which keeps constant the number of neighbors, rib, of a 
given particle. Although the simulations described below were 
obtained for rib - 100, the same calculations carried out us- 
ing ni, = 36 neighbors did not provide significant differences. 
Corrections due to the gradient of the smoothing length were ex- 
plicitly included in the IADq formulation as it comes from the 
Lagrange-Euler variational principle. 

In Table 1 we show a summary of energy and momentum 
conservation properties for the different calculated models. The 
tensor calculations are always more able to ensure momentum 
and energy conservation than the standard ones for the same 
elapsed time. The lack of energy conservation during the point- 
like Sedov explosion in all methods is due to the hard initial 
conditions. In contrast, the very good conservation of momen- 
tum in the supersonic numerical experiments is partially due to 
the spherical symmetry imposed to the initial configurations. 
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4.1. Kelvin-Helmholtz instability 

The Kelvin-Helmholtz instability is a well-known test to check 
the ability of any hydrodynamic code to handle subsonic pertur- 
bations. This instability appears when there is a sufficient veloc- 
ity shear in the interface layer between two fluids with different 
densities. Small perturbations of the velocity field in the orthog- 
onal direction to the interface grow, leading to a mixing of both 
fluids. This is usually simulated in a box with periodic boundary 
conditions where two fluid regions are defined with densities p\ 
and p2 respectively. Both layers have opposite parallel veloci- 
ties leading to a shear discontinuity in the contact interface. To 
develop the instability, a small perturbation is seeded in the in- 
terface as a sinusoidal mode of length scale A. 

In our case, we have simulated a central band of high den- 
sity fluid (pi) moving in a low-density medium (P2) in a squared 
lattice of 1 cm side in the XY plane using N=62500 particles. 
The mass of the particles was arranged to obtain the corr e ct den - 
sity profile following a ramp function [Robertson et al.l (|2010). 
In this way, we smoothed the interface density jump to make it 
comparable to the SPH resolution using 

f(y) = A i +exp 2^ 1+exp 2(o^ ' (38) 

where A is a normalization constant and Ay = 0.05 cm. 
Therefore, the density profile was given by 

p(y)=P2 + (Pi-P2)/(y), (39) 

where p\ —2 g.cuT 3 and p2 = 1 g.cm -3 . 

The seed of the perturbation is obtained using a sinusoidal 
function for the v y component of the velocity field. For the initial 
velocity, we then have 

v. x (y) -v 2 + (vi - v 2 ) f(y) , v y (x) = Av y sin (mix) , (40) 

where we assume n — 2 and Av v = 0.1 cm.s -1 . Also note 
that the v x component has been smoothed using the same ramp 
function used for the density, and that vi = 0.5 cm.s -1 and 
V2 = -0.5 cm.s 1 , which corresponds to the high and low density 
bands respectively. 

Figure [4] shows four snapshots of the growth of the Kelvin- 
Helmholtz instability at different times for the calculation using 
IADq (first row) and the standard SPH implementation (second 
row). As can be seen, the standard formulation does a poor job of 
resolving the structure of the instability. Although in the begin- 
ning the main shape is in gross agreement with the tensor simu- 
lations, the final image is blurry, the interface is not well-defined 
and the shape is incorrect. In the case of the tensor calculation, 
the instability grows cleanly and at good rate, and the definition 
of the extremes of the billows, where the finest structure appears, 
is clearly enhanced. To achieve similar results to those obtained 
with the IADq technique, different methods have been proposed 
to maintain the standard descripti on, mainly based on including 
an artificial thermal conductivity dPrice 120081) . This method pro- 
vides reliable results, but includes a new set of parameters and 
estimates, such as maximum signal velocity between two parti- 
cles to obtain a diffusion parameter. Furthermore, the inclusion 
of an artificial thermal conductivity leads to an extra dissipation 
of gradients away from discontinuities, hence it needs to design 
some means of controlling the amount of dissipation. 

In Fig. [5] we present the evolution of total linear momen- 
tum during the development of the instability. Although both 



schemes give a very satisfactory momentum conservation, we 
can see that conservation using the IADq method is at least one 
order of magnitude better than that of standard scheme. 

To test the tensor approach in a more hostile scenario, we 
diminished the value of the amplitude of the initial velocity per- 
turbation by an order of magnitude (i.e. Av v = 0.01 cm.s -1 ). 
In the third and fourth rows of Fig. [4] we show the results of 
the simulations using the configurations IADq and STD. It is 
clear that in the standard formulation the instability was unable 
to grow, while it does using IADq. Changing the geometry of 
the initial particle distribution from a square to a uniform hexag- 
onal lattice or reducing the size of the smoothing length tak- 
ing rij, = 36 neighbors did not appreciably change the results. 
We can then state, that the matrix approach has inherent virtues 
that are absent in the standard formulation, mainly related to the 
generalization of the derivation technique to a tensor expression, 
which, in some sense, diminishes the errors derived from the dis- 
cretization. 

It is also worth to mention that during the simulations car- 
ried out using the standard calculation a generalized clumping 
of particles in the low density regions often appeared. This is a 
well-known problem, called pairing or tensile-instability, where 
particles get "stuck" if r, ; - < | h owing to an unstable stress- 
strain relation that occurs typically when the second derivative 
of the kernel becomes negative and high strain is developed be- 
tween particles. This problem eventually leads to difficulties in 
resolving the distances between neighbors, incorrect interpola- 
tions, and finally large non-physical behavior of the particles 
that halts the calculation. A typical way to cope with this prob- 
lem is to adaptively c hange the kernel slop e so that it becomes 
more centrally peaked (Cabe zdn et al.l20 08). However, the IADq 
method was unaffected by this problem, hence we found that the 
tensor method seems to be less susceptible to pairing instability. 
The ability of matrix methods t o cope with pairin g instability has 
been reported by other authors (Q ger et alj20 07). and they merit 
further investigation in the context of IADq. 

4.2. Rayleigh-Taylor instability 

The simulation of the growth of the Rayleigh-Taylor instability 
in a stratified fluid in the presence of gravity is a classical test of 
subsonic fluid dynamics. In its simplest form, a box containing 
two fluids with different densities separated by a sharp transi - 
tion region is placed inside a gravitational field dYoungslll984h . 
If the denser fluid is located on top of the lighter, the system 
becomes physically unstable and any small perturbation of the 
interface leads to fluid overturn. The denser fluid sinks and the 
lighter one rises, producing characteristic structures known as 
spikes and bubbles respectively. After a linear phase of growth, 
which can be studied analytically, the system enters a complex 
non-linear regime. In real fluids, however, the viscosity of the 
fluid prevents the growth of the instability at short wavelengths, 
thus viscosity tends to dampen or even to completely suppress 
the instability. One serious problem of numerical schemes incor- 
porating the artificial viscosity formulation is that they introduce 
too much viscosity into the system, compromising the growth 
of instabilities. As we later show, the integral approach to the 
derivative, Eqs. ( 1331 ) and ( T37I ). provides a substantial improve- 
ment in the simulation of the RT phenomenon when the initial 
perturbation amplitude is small. 

A sample of N=62 500 particles was distributed in a squared 
lattice of 1 cm side. The mass of the particles was conveniently 
arranged to reproduce the density profile given by 
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P(x,y) = 



P2 
[Pi 



y > Ay 

-Ay < y < Ay 
y < -Ay 



(41) 



where y is the vertical coordinate with the origin at the interface 
and Ay is the width of the transition zone between both fluids. 
The density in the two zones were set to p\ — 1 g cirT 3 and 
P2 = 2 g cm 4 respectively, and the size of the transition zone 
was taken to be Ay = 0.05 cm. 

The integration of the hydrostatic equilibrium equation gives 
the pressure distribution along the fluid 



{-gP2(y-k)\ 



P(x,y) 



y > Ay 



-g 



1 



pi (y - -)+ 

( P2 A ~ Pl )(y 2 - Ay 2 )+ -Ay <y< Ay 



4Ay 

( — z — )(y- Ay) 



(42) 



-piy); 



y < -Ay 



where g is the gravity, here of value g=0.5 cm s~ 2 . With this set- 
ting, the speed of sound at the inter-phase is c s = 0.7 cm s _1 . A 
perturbation in the boundary layer with A = 0.25 cm was seeded 
by giving an initial small ver tical veloc ity v y to the particles ac- 
cording to the prescription (lAbel ll20TTl) 



AVy 

v y (x,y) = — <l +cos 



8;r 



1 

X+ 4 



1 + cos 



5n\y-- 



(43) 

and v y was set to zero for y positions above 0.7 cm and below 
0.3 cm. The value of velocity perturbation was set to Av y = 
0.01 cm s" 1 . In Fig.|6]we show several snapshots of the develop- 
ment of the Rayleigh Taylor instability calculated using IADq. 
As we can see, the RT instability is able to grow after a few 
tenths of a second. When the same calculation was attempted 
using the standard SPH scheme to compute the acceleration, the 
results were of much lower quality, the instability was totally 
damped and there was no development of the RT fingers. As in 
the case of the Kelvin-Helmholtz instability, changing the initial 
particle setting from the square to an hexagonal regular lattice or 
reducing the size of the smoothing length did not qualitatively 
alt er the above conclusions. A similar conclusion was reported 
bv lAbei1(l20Tl who used a variant of standard SPH, albeit not 
fully conservative, to simulate the grow of the RT instability us- 
ing a similar initial setting of the experiment. This negative result 
does not mean that standard SPH is unable to simulate this phe- 
nomena because it could handle it after a careful initial setting 
that minimizes the numerical noise. Nevertheless, the results of 
our simulations indicate that the tensor method is less dependent 
on the particular geometry of the initial lattice. 

In Fig. [7] we show the evolution of the center of mass of 
the spikes y s from the initial contact surface. According to 
the theory, the evolution is exponential during the initial stage 
£(f) = £o exp [T(t - to)], roughly until the initial perturbation £o 
at fo (here we took f = 1.1 s) has grown to a size =2 l/k 
where k is the wave number. In our tests, 0.04 cm. The the- 
oretical growth rate F, during the exponential phase is given by 
T f = y/A, kg, where A, = (P2 - Pi )/(p2 +Pi ) is the Atwood num- 
ber. The exponential phase is followed by a brief second stage in 



which the penetration of the dense fluid evolves as y s oc gt 2 un- 
til drag takes over and both falling spikes and rising bubbles 
reach a limiting speed. According to Fig. |7j all these features 
are present in the simulations. Nevertheless, there was only a 
qualitative agreement with the analytical growth rate during the 
initial exponential phase because the value of F inferred from 
the simulations was almost one half of the theoretical value. This 
is unsurprising because a reliable quantitative estimation during 
the initial stage needs both a higher resolution as well as a better 
initial setting than that we are currently using. In particular, the 
transition zone between light and dense fluids has to be sharper 
to correctly represent the Atwood number in this stage. An ana- 
lytical expression for the terminal velocity of a rising bubble in 
a tube was obtained bv lLavzerl d 19551) . v h = 0.51 yfX^i where 
/ is the radius of the tube. Taking I - A - 0.25 cm, we get 
Vb = 0.104 cm.s -1 in good agreement to the numerical value 
deduced from Fig. [7] 



4.3. The wall heating shock test 

We now consider simulations of supersonic events. In the test 
known as the wall shock problem, a spherical or cylindric su- 
personic stream of gas is launched towards its geometrical cen- 
ter forming a highly compressed region. For these geometries, 
simulations can be compared to the predictions of an analytical 
approach to the evolution of th ermodyn amic variables as a func- 
tion of the initial conditions dNohl ll987). Although the gross 
features of the event are correctly captured by SPH simulations, 
it is well-known that schemes based on artificial viscosity have 
difficulties in providing a detailed description of the wall heating 
shock test. The reason is that AV spreads the shock over several 
computational cells, which induces a non-physical increase in 
the internal energy ahead of the shock. For converging flows, a 
large artificial spike in internal energy is observed at the conver- 
gence center. As a consequence, a profound dip in the density 
profile appears to keep the pressure smooth. However, the in- 
clusion of a good amount of artificial viscosity is mandatory to 
providing a successful description of the shock. 

Our purpose behind this test was to discern whether the in- 
tegral route to calculating derivatives can provide more reliable 
results than the standard procedure. Therefore, we use IADq and 
Eqs. $15[ , d33l and (l37l l to carry out the simulations. 

A sample of N=57 600 particles of the same mass were 
evenly spread across a lattice to ensure that the density was ho- 
mogeneous, p = 1 g cm 3 . The initial pressure and internal en- 
ergy of particles were negligible. The system was imploded by 
imposing a spherically symmetric velocity field vy = -1 cm s . 
In Fig. [8] we show the density and radial velocity profiles at 
t=0.3 s, when the shock is well-developed. Both profiles follow a 
simple pattern that remains close to the analytical estimations. In 
the central region, r < 0. 1 cm, a plateau of compressed material 
forms, with a characteristic density of p s ^ 14.5 g cirT 3 a little 
lower than t he analytical value of p = 16 g crrT 3 for y — 5/3 
dNoh II 19871) . In this inner region matter remains stagnated with 
a velocity close to zero. Outside the plateau the density abruptly 
declines through the shock front trying to regain its initial value, 
whereas the velocity increases to reach v r — -1 cm s at the 
incoming, but still unshocked, matter. As we can see in Fig. [8] 
the simulation using the tensor approach is of higher quality than 
that of standard derivatives. In particular, the numerical oscilla- 
tions in the post-shock region in both, density and velocity, were 
considerably smaller in the tensor formulation. Nevertheless, the 



7 



Garcfa-Senz, Cabezon and Escartm: Improving SPH using IADq 



dip in the density profile and the maximum value achieved by 
density were similar in both calculations. 



4.4. Sedov test 

In the Sedov test, the evolution of a shock wave front born as 
a consequence of a point-like explosion is studied as it prop- 
agates in a homogeneous medium. The problem of an intense 
explosion in a gas is a standard test for hydrocodes that is of rel- 
evance to astrophysics, where it is not rare to find strong shocks 
in many scenarios involving fluid motions at high velocity. The 
theoretical solution was found by L.I. Sedov by applying self- 
simil ar methods an d dimensional analysis for different geome- 
tries (Sedov 1959). In its simplest formulation, the Sedov prob- 
lem has an initially cold gas at rest. At t = s there is a point 
explosion at the origin that in Sedov ( 1959) was treated as an in- 
stantaneous release of energy at the origin, and assumed that the 
background material through which the expanding gas sweeps 
behaves like a perfect gas. For these initial conditions there are 
precise, albeit algebraically complicated, analytical expressions 
for the fluid variables. In the case of SPH, the artificial viscos- 
ity smears the shock over 2-3 times the smoothing-length. As 
a consequence, the density jump across the shock front is al- 
ways smaller than the factor of four predicted by the theory for 
y = 5/3. In more than one dimension, resolution issues are cru- 
cial not only to resolve the peak of the blast wave but also to 
reproduce the correct post-shock variables downstream and the 
structure of the rarefied tail close to the origin. We want to in- 
vestigate the extent to which the use of the integral approach to 
the derivative can improve the results of the SPH simulations. 

The initial setting was similar to that of the wall heating test 
but this time the velocity was set to zero everywhere at t=0 s, to 
ensure that the initial configuration is in equilibrium. To create 
the explosion, a large amount of internal energy was instanta- 
neously released in a small region around the center of the box. 
To smooth the initial discontinuity, we took an initial pressure 
step that decays as a Gaussian function 



P(r) = P 2 + (P l -P 2 )exp 



~7a 



(44) 



where Pi and P 2 are the pressures in zones left and right of the 
step and cr sets the width of the pressure decay. In the present 
simulations, P\ = 10 4 dyn crrT 2 and P 2 = 1 dyn crrT 2 and 
o 2 = 36 cm 2 , which smooths the internal energy step over about 
5-6 times the smoothing length. In Fig. [9] we depict the den- 
sity and velocity profiles once the self-similar state has been 
achieved. As we can see the profiles calculated using the tensor 
and standard approaches to calculating gradients do not differ so 
much. The peak in the density profile is about 80% of what is 
expected from a strong shock moving through a y — 5/3 gas but 
the profile in the post-shock zone is smoother in the calculation 
using IADq. The velocity profile is also similar but the tensor cal- 
culation showed fewer oscillations and a more ordered behavior 
in the post-shock tail. These results agree with our main con- 
clusions given in the previous test dealing with the wall heating 
shock, reinforcing the idea that for identical initial settings IADq 
provides more reliable results than standard derivatives also for 
supersonic events. 



5. Handling linear phenomena: momentum and 
energy equations using IAD 

Ensuring both exact linear interpolation and perfect momentum 
and energy conservation using either full IAD or other similar 
schemes, seems to be difficult. In this respect, the best physical 
systems to try to conciliate both things are pure linear systems 
for which linear acoustic wave propagation is a representative 
example. A suitable alternative SPH formulation can be found 
starting from the differential acceleration equation 



1, 
P 



-VP 



(45) 



which, using Eq. (fTUl i. can be easily accommodated within the 
SPH formalism 



b \Pa P b ) 



ab ■ 



(46) 



This equation can be deduced independently using the vari- 
ational principle and, provided the kernel is symmetrized W a b = 
Q.5(W a b{h a ) + W ab (h b )), it conserves linear and angular momen- 
tum exactly. If instead of Eq. dTOb we use Eq. ^ to compute Eq. 
( l45l l. the resulting equation is 



V|£) + £vp 



Pi 



Pi) 



VW ab + 2 P a V —VW ab . 

Pa Pb 



(47) 



This particular form of the acceleration equation has several 
interesting features: 1) In the continuum limit it reduces to the 
standard expression, Eq. d46b . because the last term in the equa- 
tion vanishes. 2) As shown in Sect. 3.2, an isobaric system with 
smooth density gradients remains in equilibrium with negligible 
acceleration for a long time (see Fig. [3]). 

However, it is also evident that the inclusion of the last term 
in Eq. ( l47l i breaks its symmetry. Still, in the linear limit, momen- 
tum is conserved to a high extent as can be easily demonstrated 
in ID. We first expand the magnitudes P/p 2 in Eq. (l47l i assuming 
a smooth density gradient 

P a P a ( 1 _*P\ P± a p b L , Ap 



, , , - ,1 + — , (48) 

Pa PaPh \ Pb) P% PaPb \ Pa j 

where Ap = p a —p b and |Ap|/p << 1. Inserting Eq. ( l48b into Eq. 
( l47l i and neglecting higher order terms, the acceleration is given 
by 



Zb%(^)(Xb-Xa)W a b 



Zb - Xa) 2 W< 



(49) 



ah 



If we assume a linear pressure profile within the kernel domain 
of the particle 



(8P\ 

P b = P a + y^J (x b - x a ) , 



the acceleration becomes 



1 tdP 



Pa \ dx ) a 



(50) 



(51) 
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which is the one-dimensional Newton equation, thus leading 
to an exact value for the acceleration. In consequence, global 
momentum is well-preserved provided that the system remains 
in the linear regime and the density gradients are smooth. 
According to this discussion, we complete the acceleration equa- 
tion obtained in Sect, fusing the IADq scheme with a corrective 
term 



V< { Pa Pb , 

~ yPa \ 

flfli 3\.i,ab ^i,ab(h a ) , 

PaPb J 



(52) 



where < *P < 2 is a parameter controlling the strength of the 
applied correction. For *P = 0, the already checked IADq scheme 
results, whereas for T = 2we have full linear IAD. 

To take into account the corrective term in the energy equa- 
tion we include the power released/absorbed by such term in Eq. 

G3 



ldu\ _ y V mi ( v _ v ) . j Jli,ab(ha) + 



2 PaPb 



(53) 



such that total energy is exactly conserved. For very subsonic 
movements, the energetic contribution of the correction term is 
small. For = 0, Eqs. d52l and d53T l reduce to the fully conser- 
vative IADq scheme. Alternatively, the case of *F = 2 could be 
used to handle linear systems with smooth density gradients as 
it shows enhanced interpolation abilities and good conservative 
properties. 

Sound wave propagation is an example of a system that 
evolves in the linear limit of Euler equations. It is therefore a 
test ideally suited to highlighting the potential advantages of the 
IAD scheme. A homogeneous system with po = 1 g crrr 3 was 
simulated using N=62500 particles uniformly distributed in a 
lattice of size 1 cm. The initial values of both pressure and den- 
sity were set to one. A sample of n p = 200 particles inside a 
circle of radius R p - 0.032 cm located at the center of mass 
of the lattice was obliged to oscillate with small amplitude and 
period P. This compact set of particles emulated the effect of 
an external oscillating piston moving adiabatically onto the ini- 
tially unperturbed gas. As the piston moves, it launches regu- 
lar trains of circular waves that propagate at the sound speed 
c s - ^JyPjp - 1-29 cm s _1 for y = 5/3. We want to check the 
ability of the different aforementioned schemes to handle with 
this problem. 

Particles belonging to the piston move homologously, fol- 
lowing the harmonic oscillator law 



v r (r, i) = v™(r) cos(wf) , 



(54) 



where d stands for the angular frequency, and the maximum 
value of the radial velocity of a particle located at distance r 
from the center is 



v?(r) = V?{R p ) 



Rr 



(55) 



was set to P=0.05 s so that around eight complete waves were 
launched before the first one reached the limits of the system. 
The value of the maximum velocity at the piston head was set to 
v''<(R p ) = 0.16 cms- 1 . 

The results of the simulations for the three SPH approaches: 
IAD, IADq and Standard are summarized in Fig. [TO] which rep- 
resents the radial velocity profile of the gas at time t=0.35 s. It is 
quite evident that the spherical symmetry was better preserved 
during the IAD calculation, while there is a progressive degra- 
dation of the symmetry in the IADq and standard calculations 
respectively. The worst case corresponds to the standard scheme 
but even in this case spherical symmetry is reasonably preserved. 
A more quantitative analysis can be done making use of the an- 
alytical solution for diverging circular waves, usually known as 
waves in a membrane in the literature. The solution for free har- 
monic traveling waves in a circular membrane is written as 



A(r, t) = A,„J(u)cos(wt) , 



(56) 



The resolution places limits on both the period and the max- 
imum velocity at the piston head v'"(R p ). The value of the period 



where u = kr, k is the wave number and J(u) is the Bessel func- 
tion with radial and angular modes m=n=0. The Bessel func- 
tions are constrained by 7(0) = 1 at the center of the wave 
train. The profile for the radial velocity given by Eq. d56b with 
A m = v'"(R p ) and the comparison with the numerical profiles 
obtained using the different schemes are also shown in Fig. \W\ 
The analytical solution reminds one of a damped cosine wave 
propagating at the sound speed c s = 1.29 cm. . On the whole, 
the three schemes were able to correctly track the evolution of 
the waves. However, a close inspection reveals small differences 
among them. The best value of wave velocity corresponds to 
the IADq scheme (c'" d ° = 1.285 cm.s -1 ) followed by the stan- 
dard (cf d = 1.267 cm.s -1 ), whereas full IAD falls a bit short 
(Cj = 1.25 cm.s -1 ). Even though the damping in the three 
schemes is always larger than that of the analytical solution, 
it seems that the worst case also corresponds to the fully IAD 
scheme built setting *P = 2 in Eqs. ( 1521 and ( 1531 . 

The simulation of linearized systems represents the most fa- 
vored situation to apply exact linear interpolation schemes such 
as IAD because, as shown above, conservation of momentum 
and energy is probably as good as in the standard formulation. 
Nonetheless, the numerical experiments with the acoustic waves 
suggested limited advantages of using these schemes. In fact, 
some details of the wave evolution are better represented by 
IADq or even by the standard scheme. Only the spherical sym- 
metry was more accurately preserved in the experiments using 
IAD. 



6. Applying IADq to astrophysics: the stability of a 
Sun-like polytrope 

Applying Eqs. ([TBI . ( 1331 and (l3~7l to astrophysics is straight- 
forward by just adding the gravitational acceleration to the mo- 
mentum equation. W e approached gravity us ing a multipolar ex- 
pansion of the force (Hernqu ist & Katzll 19891) . up to quadrupole 
contributions. The 3D structure of a Sun-like polytrope was sim- 
ulated using a sample of N = 10 5 particles with equal mass. 
To achieve the equilibrium, we proceeded in three stages. First, 
the radial coordinate of each particle was set following the ID 
density profile of the polytrope with its angular position deter- 
mined at random. In the second stage, we allowed the particle 
sample to relax under the action of pressure and gravity forces 
but the movement of the particles was constrained to keep their 
distance to the center constant. In a third stage, we allowed the 
particle sample free to approach the equilibrium configuration. 
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We followed the third stage using both IADq and STD schemes 
with the AV coefficients set to a = 1, J3 = 2. The main goal of 
this test is to show that IADq (and probably other matrix meth- 
ods) can be used to simulate astrophysical scenarios with good 
results, excellent conservation properties and with low computa- 
tional penalty. 

In Fig. QT| (upper row) we depict the evolution of density for 
two particles initially located at the center and the surface of the 
star respectively. As we can see, these particles display the typ- 
ical oscillatory behavior in both regions. As expected, the am- 
plitude of the oscillations is small close to the center and larger 
at the surface. The relaxation towards equilibrium is slower for 
the tensor method but the structure of the star after t = 2 h of 
relaxation is very similar. With the exception of a tiny region 
close to the surface, we verified that the gradient of pressure 
matched gravity along the poly trope in both calculations. In Fig. 
[Tn (lower row) we show the evolution of the kinetic energy and 
the instantaneous position of the center of mass during the re- 
laxation process. The evolution is similar in both calculations 
but the level of kinetic energy is always a bit higher when iADo 
is used. The energy conservation after two hours of evolution 
is rather good for both schemes =* 3 10~ 5 . In the bottom right 
of Fig.[TT]we show the evolution of the moduli of the center of 
mass position. We see that momentum is not perfectly conserved 
during the simulation in both simulations, mainly owing to the 
multipolar approach to the gravitational force. Nevertheless, the 
tensor scheme provides a better conservation of momentum than 
the standard SPH calculation. 

The higher level of kinetic energy and the slower relaxation 
rate towards complete mechanical equilibrium shown by iADo 
suggest that for disordered systems the numerical noise could be 
higher for the tensor method. As shown in Fig. QT| the recalcu- 
lation of the evolution of the polytrope using larger values of the 
AV coefficients, a - \.5,f3 — 3 led to a significant reduction in 
the level of spurious kinetic energy, while the evolution of the 
density at the center and surface of the star remained unaltered. 
This implies that the numerical noise is the main agent respon- 
sible for the excess of kinetic energy seen during the relaxation 
process. 

To explore the performance of the algorithm, the tolerance 
parameter 8, controlling the accuracy of the multipolar expan- 
sion, was varied from < 9 < 1, where 8 = corresponds to 
direct particle to particle interaction. This benchmarking anal- 
ysis is summarized in Fig. Q~2] We see that for standard val- 
ues of the tolerance parameter, 9 ^ 0.7, widely used in current 
SPH calculations, the computational overload is around 40%. 
This overload rapidly decays as 6 — » 0, as expected. However, 
the number of operatio ns in the multipolar calcu lation of grav- 
ity scales as NlogN dHernquist & Katj [T989). We therefore 
expect a reduction in the computational overload with increas- 
ing numbers of particles. Nevertheless, we note that all these 
calculations were performed using a serial code. An invaluable 
property of SPH is that only a single tree-walk is needed to find 
the neighbors (and calculate gravity if needed). Once this infor- 
mation is available the remaining calculations (density, momen- 
tum and energy equations, IADq terms) can be performed with 
linked-lists that can be directly parallelized. Following that idea, 
we found that for a 3D simulation of 2 ■ 10 5 particles in a stan- 
dard 4-core desktop computer, the gravity calculation (which re- 
mains serial) takes around 30 times more wall-clock time than 
the rest of the calculations together (parallelized with OMP us- 
ing 4 cores). Knowing this, the computational overload due to 
the IADq calculation remains sub-dominant. 



7. Conclusions 

We have presented and verified a novel procedure (IAD) to eval- 
uate gradients in the context of SPH simulations. The mecha- 
nism for calculating derivatives relies on an integral approach, 
which ensures that the derivative of a linear function is ex- 
actly obtained, even after transforming integrals to summations. 
The drawback is the greater algebraic complexity because gradi- 
ents are estimated using the tensor expression given by Eq. 
Nevertheless, our new scheme is not significantly more complex 
than the standard one, easy to implement and includes the clas- 
sical formulation as a particular case. A shortcoming of our new 
method is that it leads to non-conservative movement equations. 
Thus, we have developed and tested a restricted interpretation 
of the method, referred to as IADq, which sacrifices exact lin- 
ear interpolation to ensure a perfect conservation of momentum 
and energy. The analytical considerations and numerical exper- 
iments described in Sects. [2] and |4] strongly indicate that IADo 
behaves better than the standard method in computing gradi- 
ents. The modified momentum equation obeying these deriva- 
tive rules was developed in Sect. [3] resulting in Eq. ( l29b . Since 
the movement equation was obtained from the Euler-Lagrange 
variational principle assuming isentropic evolution, the ensuing 
equation was fully conservative and explicitly included the Vh 
terms. A conservative energy equation compatible with the mo- 
mentum equation was also developed. To handle with shocks 
the scheme was completed by including of the standard artifi- 
cial viscosity formalism resulting in Eqs. ( f33b and (l37l i. which 
in addition to the density equation in Eq. d!51 l. summarizes the 
mathematical formalism linked to iADo. 

The formulation of SPH using matrix methods based on the 
variational approach dBonet & Lokll999l) has been used in CFD 
to successfully simulate a variety of problems (generally in two 
dimensions) from fluids to the impact and fracture of solid bod- 
ies. Nevertheless, none of these schemes are able to simulta- 
neously achieve re-normalization, exact momentum and energy 
conservation, perfect linear interpolation and implici tly include 
the gr adient of the smoothing length in the equations (lOger et al.l 
2007). Nevertheless, iADo is probably the optimal formulation 
because it fulfills almost all the above requirements with a mod- 
erate computational overload. 

Four tests were performed in two dimensions to verify the 
performance of the method (see Table 1), two of them in connec- 
tion to bi-dimensional subsonic hydrodynamic evolution (KH 
and RT instabilities) and the other two related to the descrip- 
tion of highly supersonic phenomena such as the wall heating 
shock and the Sedov tests. In all cases the performance of the 
new scheme was superior, although in the supersonic numeri- 
cal experiments the improvement was modest. The tests of the 
growth of the Kelvin-Helmholtz and Rayleigh-Taylor instabil- 
ities clearly showed the power of the IADo scheme because 
no growth at all was seen when the standard scheme was used 
with small initial perturbations, whereas instabilities were able 
to grow when iADo was used. This does not of course mean that 
the standard SPH technique cannot give a satisfactory answer to 
these problems but simply states that for the same initial condi- 
tions the tensor method seems to be more stable and sensitive to 
small perturbations. Therefore, these results indicate that simu- 
lations using the IADq scheme are less dependent on the initial 
setting of the particles. As expected of a tensor method, this ad- 
vantage could probably be reinforced as the dimensionality in- 
creases. 

In Sect. [5] we devised and discussed a slightly different for- 
mulation, which combines exact linear interpolation and good 
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conservative properties, given by Eqs. ( f52b and d53l l. The main 
point here is that the fully IAD method would be only useful 
for describing linearized systems with smooth density gradients. 
Nevertheless, a single parameter allows us to easily switch from 
IAD to IADq and in this sense the formulation given in Sect.[5]is 
more general. As sound waves are the prototype of a linearized 
system we have used the new scheme to simulate the 2D propa- 
gation of acoustic waves. The simulations indicate that there are 
no particular advantages of using IAD instead of either IADq or 
standard SPH, because only the symmetry was better preserved 
during the evolution of the wave trains. 

A preliminary application of IADq to astrophysics was dis- 
cussed in Sect. [6] in connection to the stability of a Sun-like star 
described by a poly tropic EOS. The main goal was to demon- 
strate that matrix SPH methods can be used to handle 3D self- 
gravitating bodies with satisfactory results and to explore the real 
computational overload when long-range forces need to be com- 
puted. For a reasonable accuracy in the calculation of gravity, 
we have estimated a computational overload < 50% for a serial 
code, with respect to that of standard SPH, when a multipolar ex- 
pansion is used to calculate the gravitational force. The overload 
becomes negligible for a precise particle-to-particle calculation 
of gravity scaling as N 2 12. 

To summarize we could say that simulations carried out us- 
ing the SPH scheme obtained with the IADq approach to the gra- 
dients always led to improved results with respect to standard 
SPH. As the new scheme is both fully conservative and more 
precise in making interpolations, it could be an alternative to the 
standard technique in handling systems subjected to small per- 
turbations. This conclusion is supported by the results of our nu- 
merical study of the growth of Kelvin-Helmholtz and Rayleigh- 
Taylor instabilities. In addition, our simulations of supersonic 
phenomena also improved when the tensor approach was used. 
The main drawback of this method is that it increases the com- 
putational overload, but one has to keep in mind that there is 
not always a linear relationship between algebraic complexity 
and computational charge. This could be true for hydrocodes 
that incorporate time-consuming physics or when the SPH al- 
gorithm built with IADq allows longer time steps. Even more, 
using linked-lists a direct parallelization of the method is possi- 
ble as the calculations are carried out of the tree-walk, keeping 
the computational overload of the IADo almost negligible com- 
pared to more time-consuming sections of the code. 

Although the presented results are encouraging more work 
needs to be done to confirm and extend the conclusions of our 
proposal. For the most part, the tests presented in this work to 
validate the scheme were carried out in 2D boxes using well- 
ordered initial models. The simulation of realistic astrophysical 
scenarios generally involves, however, a quite different numeri- 
cal setting. Many of these calculations have to be performed in 
3D using a random-like initial particle distribution and incorpo- 
rating the long-range gravitational force. A detailed compara- 
tive analysis of the ability of IADq and standard SPH to cope 
with these scenarios is beyond the scope of the present work and 
is left to a forthcoming publication. Nonetheless, we can sug- 
gest a couple of areas of difficulty in the tensor formulation that 
will probably come up in astrophysical 3D applications of the 
method: 1) As suggested in Sect. [6] the tensor method might 
be more sensitive to the disorder of the particles and display 
a higher level of numerical noise than the standard scheme. 2) 
Free boundary conditions could be more difficult to handle using 
IADq because the simplification hypothesis assumed in Eq. (fT4l > 
does not hold at the edge of the system. The difficulty in sim- 
ulating sharp boundaries with matrix methods is a well-known 



problem of MLS schemes. According to lOger et all (120071) . it 
can be solved by taking special conditions at the limits of the 
system. This is not, however, a big concern in astrophysics be- 
cause astronomical bodies never have abrupt boundaries. Thus, 
the stability test of the Sun-like star discussed in Sect.[6]does not 
reveal that the particles located close to the surface have peculiar 
behaviors. Although the initial imbalance between gravity and 
pressure gradient is more pronounced than in the standard for- 
malism, the ensuing evolution towards equilibrium is not much 
different. 
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Fig. 1. Relative error e in the first derivative of density p(x) = 
(1 + x) g cm -3 as a function of the smoothing length (in 
interparticle units) calculated using both tensor and standard 
SPH schemes. Crosses are for direct EilEs. derivative estimation. 
Upper left is for the error when the Gaussian kernel is used. 
The upper right picture shows the contribution of several error 
sources to the total error e. Bottom pictures are the same but for 
the cubic-spline kernel. 
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Fig. 2. Relative error e in the first derivative of linear functions 
calculated in 2D as a function of the smoothing length. Upper 
rows depict the value of e for p(x,y) = (1 + x) g cirT 3 ob- 
tained using the Gaussian (left) and the cubic spline (right) ker- 
nels. Bottom rows are the same but for the function p(x,y) = 
1 + yjx 2 + y 2 g cm 4 . Error profiles with crosses are for the di- 
rect derivative estimation. 
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Fig. 3. Evolution of magnitude <t p ^ - ^J ^ P '' N of an inho- 
mogeneous 2D system, initially in hydrostatic equilibrium, cal- 
culated using the different SPH schemes mentioned in the text. 




Fig. 4. Evolution of the Kelvin-Helmholtz instability. The snap- 
shots show the color map of density at times t = 0.1, 1.0, 2.0 
and 3.0 s for methods IADo (first row) and standard SPH (second 
row), with Av y = 0.1 cm s -1 . Times t = 1.0, 3.0, 4.0 and 5.0 s 
for methods IADq (third row) and standard SPH (fourth row), 
with Av\. = 0.01 cm s _1 . 
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Fig. 5. Evolution of the absolute deviation from initial linear mo- 
mentum during the development of the Kelvin-Helmholtz insta- 
bility. Figure on the left is for the x component of total mo- 
mentum calculated using IADq and standard (STD) schemes, 
whereas figure on the right is the same for the y component. 
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Fig. 6. Development of the Rayleigh-Taylor instability for 
Atwood number 1/3 calculated using the IADq scheme. The 
snapshots show the color map of density at times t = 
0.4, 4.2, 5.1 and 5.7 s. 
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Fig. 7. Evolution of the center of mass of the spikes during the 
growth of the Rayleigh-Taylor instability depicted in Fig. [6] The 
velocity of the center of mass of the bubbles, Vb, and spikes, v s , 
is also shown. For t > 4 s, a limiting speed is reached. 
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Fig. 9. Density and velocity profiles during the Sedov test. 
Figure on the left is for IADq calculation, whereas the figure on 
the right is for the simulation using the standard SPH scheme. 
Density is normalized to 4 g cm 4 . 



Fig. 10. Acoustic wave profiles at t=0.35 s calculated using 
IAD, IADq and standard SPH (STD). The waves were gener- 
ated by the periodic displacement of a circular piston of size 
0.035 cm located at the center of the lattice. The continuum 
line is the analytical solution calculated taking a sound speed 
of c s = 1.29 cm s _1 . 
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Fig. 11. Evolution towards stability of a Sun-like star approached 
by a polytrope. Calculations were carried out using JADo and the 
standard, STD, schemes. Labels IADq (1 ) and IADo (2) refer to 
calculations with the AV parameters set to a = 1, (3 — 2 and 
a — 1.5, p — 3 respectively. 
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Fig. 8. Density and velocity profiles during the wall heating 
shock test. Figure on the left is for IADq calculation, whereas 
figure on the right is for the simulation using the standard SPH 
scheme. Density is normalized to 16 g cm" 3 . 
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Fig. 12. Benchmarking of IADq versus STD. Parameter 6 is the 
tolerance parameter assumed in the multipolar calculation of 
gravity. 
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Table 1. Conservation properties of the computed models: Kelvin-Helmholtz (KH), Rayleigh-Taylor (RT), the wall heating shock 
(WH) and Sedov (SED). Conservation of momentum is tracked by normalizing the deviation of the center of mass over the charac- 
teristic size of the system R. See Fig.[5]for the temporal evolution of the momentum in the KH case. The method used to calculate 
the derivatives IADq or standard, STD, is indicated. 
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